Theoretical Design of a Shallow Donor in Diamond by Lithium-Nitrogen Codoping 
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We propose a new substitutional impurity complex in diamond composed of a lithium atom that 
is tetrahedrally coordinated by four nitrogen atoms (LiN4). Density functional calculations are 
consistent with the hydrogenic impurity model, both supporting the prediction that this complex 
is a shallow donor with an activation energy of 0.27 ± 0.06 eV. Three paths to the experimental 
realization of the LiN4 complex in diamond are proposed and theoretically analyzed. 

PACS numbers: 71.55.Cn,71.15.Mb,61.72.Bb 



C/5 



-a 
c 
o 

(N 
> 

in 

(N 

in 
o 

(N 



X 



With respect to multiple figures of merit that estimate 
semiconductor performance in high-power electronics [l[ , 
diamond is the best of all known semiconductors. Steady 
progress is being made to realize diamond electronics by 
improving the quality and reducing the cost of single- 
crystal diamond films made by chemical vapor deposition 
(CVD) [2]. Substitutional boron doping has succeeded in 
producing p-type diamond and can be incorporated up 
to concentrations sufficient to cross the metal-insulator 
transition and produce superconductivity [3j. However, 
no donor impurity has been incorporated into single- 
crystal diamond with sufficiently small activation energy 
and high concentration to produce an n-type semicon- 
ductor suitable for high-power applications [4|. 

In this letter, we propose a new substitutional donor 
complex in diamond composed of lithium tetrahedrally 
coordinated by nitrogen and report on a theoretical study 
of its activation and formation. Favorable properties of 
LiN4 can be inferred from properties of similar structures. 
Lithium tetraamine [5j is a metal composed of Li(NH3)4 
molecules that are locally isostructural and isoelectronic 
to LAN4 in diamond. Each Li(NH3)4 molecule donates 
an electron to a metallic state permeating the intersti- 
tial region between molecules. If this effect persists for 
dilute UN4 in diamond, it should produce shallow donor 
states with a small activation energy. A small formation 
energy is expected for LiN4 based on the high stability 
of the B center in diamond @|, which is a vacancy (V) 
tetrahedrally coordinated by nitrogen, FN4. 

Our proposal naturally follows from previous codop- 
ing proposals 0, H[ of multi-impurity complexes designed 
to prevent a carbon-nitrogen bond from breaking near a 
substitutional nitrogen impurity. The broken bond forms 
a deep, localized donor state. If the bond is preserved, the 
donor state is predicted to be shallower and more delocal- 
ized. The originally proposed BN2 and newly proposed 
LJN4 can be connected to nitrogen through a sequence 
of IN n donor complexes, CN -> BN 2 -> BeN 3 -> LiN 4 
(Fig. QJ, by reducing the valence of the central atom and 
electronically compensating with neighboring nitrogens. 



The study of other defects provides a useful reference, 
for comparison to A"N„ and as a theoretical benchmark. 
We consider two well-known substitutional impurities, 
phosphorus and boron, and two artificial defects, and 
Cf . Cf is a donor formed by adding an electron and com- 
pensating the charge by replacing a carbon and its four 
neighbors with fictional nuclei of nuclear charge Z = 6.2. 
Cf is a similarly constructed acceptor with Z = 5.8. This 
distribution of nuclear charge preserves the diamond lat- 
tice with minimal distortion and enables the formation 
of shallow defect levels. 

The minimum activation energy of a point defect is 
estimated by the hydrogenic impurity model [9(. In this 
model, activation is independent of microscopic details 
of a defect. It depends only on bulk material properties: 
the dielectric constant e and the ratio between a charge 
carrier's effective mass to* [l(| and the bare electron mass 
to. The defect loses a carrier to the nearby band edge, 
where it is Coulombically bound by 13.6— e -2 eV to the 
ionized defect. This produces a donor level below the 
conduction band edge E c or an acceptor level above the 
valence band edge E v offset by the binding energy. Using 

the model predicts defect levels 
L 0.45 eV in diamond compared 
, + 0.052 eV in silicon. 



experimental values 11 1 
at E c - 0.20 eV and E v 
to E r - 0.025 eV and E. 




FIG. 1. (color online) Predicted structures of the JfN n 
donors. First and second neighbors from X are displayed. 
Bonds are omitted for unbonded carbon-nitrogen neighbors, 
all of which are separated by 2.0 A. The remaining C-N 
bonds vary in length from 1.46-1.54 A. The X-(C/N) bonds 
lengthen from right to left on the periodic table: 1.50 — 1.60 
A for B, 1.57 - 1.65 A for Be, and 1.72 A for Li. 
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Despite its simplicity, the hydrogenic impurity model 
is empirically successful. With acceptor levels measured 
at E v + 0.37 eV in diamond and E v + 0.044 eV in silicon, 
substitutional boron accurately fits the model. The sub- 
stitutional phosphorus donor, at E c — 0.61 cV in diamond 
12| and E c — 0.045 eV in silicon, is considered a shallow 
donor in silicon but not in diamond. Shallow donor levels 
at E c — 0.23 eV have been reported in heavily deuterated 
samples of boron-doped diamond [3] , but with a lifetime 
too short for applications. 

A theoretical determination of whether LiN4 is indeed 
a shallow donor requires treatment of microscopic details. 
We use a recently proposed method [IH that decomposes 
the donor activation energy A d into a vertical ionization 
energy and structural relaxation energy, 

A D = E c + £;+ t (R£) - E° ot (R° D ) 

= [E c E D (K+)] + K t (R+) - E° ot (R D )] . (1) 

A io„i Z o Agta * 

E^t is the total energy of the system with net charge 
Q and atomic coordinates R. Ed (R) is the donor energy 
level, equal to E® ot (R) — E^ ot (R) . The equilibrium atomic 
coordinates of the donor- containing structure with net 
charge Q is denoted by R§. The corresponding expres- 
sion for an acceptor activation energy Aa is 



A A = [E a (Ra) ~ K] + [eLC&a) ~ eUK) 



(2) 



A ,o„, i 



This decomposition into elementary excitation processes 
15j , but it enables separate calculations of 

" ,1Ilze with more 



is not unique 

A rolax with total energy methods and A 



sophisticated and accurate quasiparticle methods. 

Total energies and equilibrium crystal structures are 
calculated with density functional theory (DFT) (l6. 17| 
using the Perdew-Burke-Ernzerhof (PBE) functional jl8| . 
C B and C^ 1 are modeled with alchemical pseudopotentials 
[19l | . Isolated defects are approximated with a periodic 
array of defects in a supercell of diamond. We use a 
6x6x6 face-centered cubic supercell (432 carbon atoms 
when defect-free), consistent with previous studies [20| . 
Supercells with nonzero net charge are simulated with a 
neutralizing jellium charge distribution. 

All defect structures are relaxed in their neutral and 
ionized states from multiple random perturbations of the 
ideal diamond lattice. Stable neutral XN n structures are 
shown in Fig. Q] BN2 and BeN3 also have metastable 
structures. Metastable BN2 is similar to its originally 
predicted structure Q, with two elongated C-N bonds 
of length 1.80 A rather than a single fully broken C-N 
bond. The remaining structures, including metastable 
BeN3 and all ionized defects, produce minor distortions 
in the diamond lattice that are well approximated by one 
bond length for each bonded pair of elements. 



Quasiparticle methods, unlike DFT, are constructed to 
directly model charge excitation energies such as A lonlzc . 
The state-of-the-art is the GW method 2l| , which is too 
expensive to apply to large supercells at present. Instead, 
we use the recently proposed "PBE-e" method Q , which 
approximates the quasiparticle self-energy as 

£p BE - e (r, r') = [(1 - e-> PBE (r) + « c PBE (r)] S(r r') 
- £ ->,r')y(r-r'), (3) 

with the PBE exchange u PBE (r) and correlation v PBE (r) 
potentials and a screened Fock exchange composed of the 
1-particle density matrix p(r, r') and the Coulomb kernel, 
V(r — r') = e 2 /|r — r'|. With the dielectric constant e 
set to the experimental value llj, this method produces 



a 5.48 eV band gap for diamond, comparing well to the 
experimental value of 5.5 eV. 

PBE-e is an adequate quasiparticle method for shallow 
impurity calculations because it approximates the basic 
physics of a charge carrier bound to an ionized defect. As 
in the hydrogenic impurity model, a donor state ipo(j) 
should see an effective Hartree potential originating from 
a screened ionized donor of net charge e _1 . Upon adding 
a neutral donor to pristine diamond, the Hartree poten- 
tial is modified by contributions from an updated nuclear 
charge, 5pi on (r), and an added donor electron charge, 

Sv H (r) = [v(r- r') [\^ D (r')\ 2 - 5p ion (r')} dr'. (4) 



The donor state also sees a modified effective potential 
from its self interaction in the screened Fock exchange, 



e'^D^WWir ~ r')VD(r')rfr' = 5v sX (r)^ D (v) 
6v sX (r) = -e- 1 ( V{v - r')|V^(r')| 2 dr'. (5) 



The total donor-induced potential is 8vh(v) + Sv s x{r), 
which corresponds to the bare ionized donor <5/?i on (r) and 
an effective screening cloud (e _1 — 1) |^z>(r) | 2 . A similar 
argument applies to acceptors. The screening cloud has 
the right net charge but the wrong length scale: the Bohr 
radius of the donor state is 6.3 A, but the screening length 
in diamond is estimated to be 1.5 A [13. We find PBE-e 
to be a good compromise between costly GW corrections 
to the screening length and PBE without Fock exchange, 
which has an effective screening cloud of unit charge that 
suppresses the long-range electron-impurity interaction 
and produces a donor impurity band nearly degenerate 
with the conduction band edge [17|. 

The periodic array of defects broadens defect levels into 
bands up to 0.4 eV in width. Modeling or extrapolation 
is necessary to extract an accurate activation energy. We 
use a tight-binding ansatz and a range of supercells from 
5x5x5 to 8x8x8 for extrapolation, which is described 
in detail in the Supplemental Material 17 1. 
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TABLE I. Donor and acceptor activation energies calculated 
with both the marker method [20[] and PBE-e quasiparticles, 
compared to experiment. PBE-e results are separated into 
relaxation and ionization contributions, A PBE_E = A lonlzo + 
A rclax , as in Eqs. ffl and (2j. ,5A ionizc is the RMS variance 
of extrapolation [lTj. All energies are in units of eV. 
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0.48 


0.27 
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0.00 
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0.01 
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0.27 
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1.30 


1.19 
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0.03 


0.69 
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1.67 


1.71 


0.86 


0.04 


0.85 
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0.31 


0.30 


0.30 


0.01 


0.00 
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0.37 


0.37^. 


0.31 


0.31 


0.03 


0.00 



a Mctastablc structure. 
b Experimental marker. 



Theoretical activation energies are listed in Table U 
alongside known experimental values. The PBE-e quasi- 
particle approach is compared to the semi-empirical 
marker method [20] , which calculates activation energies 
relative to an experimental "marker" impurity using PBE 
total energy differences. The marker method predicts 
larger activation energies than PBE-e. These deviations 
grow with decreasing activation energy and become as 
large as the value we are attempting to predict. This can 
be explained by derealization errors in PBE that are 
reduced in PBE-e with the addition of Fock exchange 
[lij ]. Therefore, PBE-e should be more reliable than the 
marker method as a predictor of activation energies over 
a wider energy range. Doubling the extrapolation vari- 
ance provides a wide enough confidence interval for the 
PBE-e predictions to be consistent with all experiments. 
LiN4 is shallower than the artificial shallow donor Cf 
and an activation energy of 0.27 ± 0.06 eV is consistent 
with the hydrogenic impurity model. We conclude that 
LiN4 is a shallow donor. 

Having confirmed the viability of LiN4 as a shallow 
donor in diamond, we now consider three synthesis paths. 
The first path is the diffusion of lithium into diamond 
[24| with a high concentration of B centers (VN4). The 
second path is high-pressure, high-temperature (HPHT) 
diamond synthesis [25| in the presence of lithium and 
nitrogen. The third path is CVD diamond synthesis with 
the LiN4 impurity pre formed in a seed material [26| or 
deposited molecule [271 ] . 

Nitrogen incorporates substitutional^ into diamond, 
found as an isolated center, a dimer, or clustered around 
a vacancy, VN m Q . High temperature treatment causes 
vacancies to become mobile and cluster with nitrogen to 



form mobile V"N m complexes, 
diffusion in N-free diamond \ 
(Lij) to be a mobile donor that is strongly trapped by 
vacancies. The natural extension to N-rich diamond is 
a general trapping process, Lij + VN m — > LiN m , which 
we calculate to bind at 6.88, 7.24, 8.04, 8.37, and 6.08 
eV for m = 0, . . . , 4. All sites trap strongly, but VN4 
is preferred least by Lij. The LiN m defect sequence has 
a regular trend of activity from triple acceptor (m = 0) 
to single donor (m = 4). Assuming all vacancies will be 
filled with lithium and the only acceptors are Li, LiN, and 
LiN2, then the defect concentrations n(X) must satisfy 
the inequality 



n(LiN 4 ) > 3ra(Li) + 2n(LiN) + n(LiN 2 ) 



(6) 



to prevent all LiN4 from being passivated. Therefore, 
lithium diffusion into a diamond sample is only likely 
to succeed in producing active LiN4 if the average num- 
ber of nitrogens around each vacancy in the pre-lithiatcd 
sample is greater than 3. 

HPHT synthesis of LiN4 at a detectable concentration 
requires sufficient thermodynamic stability of the com- 
plex at an accessible pressure and temperature. At zero 
temperature, we have found two pairwise decomposition 
processes that passivate shallow donor activity, 

2LiN 4 (LiN 4 ) 2 (7a) 
2LiN 4 + V — > LiN 3 + LiV^N 5 . (7b) 

The first reaction produces a LiN4 dimer with neighbor- 
ing nitrogens that break the N-N bond, which only lowers 
enthalpy below 530 GPa. The second reaction exchanges 
a nitrogen and binds an additional vacancy to the N-rich 
complex, which produces an octahedrally coordinated Li 
surrounded by CN5 . Assuming a zero chemical potential 
for V, this process lowers enthalpy at all tested pressures 
(up to 700 GPa) and has a minimum enthalpy reduction 
of 2.47 eV at 210 GPa. As a result of the process in 
Eq. {7b|, it is unlikely that LLN4 can be synthesized in 
HPHT or any other conditions that enable LiN4 and V 
to become mobile and interact with each other. 

Formation of the LiN4 complex in a CVD process from 
separate lithium and nitrogen sources is likely to be a rare 
event because it involves a coincidence of five atoms, each 
with a presumably low concentration. This problem can 
be avoided by preforming the complex within a precursor 
molecule. A suitable LiN4 precursor should be small to 
enhance volatility and simplify synthesis, closely conform 
to the diamond lattice it is to be incorporated into, and 
exist as a well-defined lithium-free molecule that strongly 
binds a lithium atom or ion. Diamondoids [29} satisfy the 
second constraint and many chelants [3(| satisfy the third 
constraint, but we propose a new analog of cyclododecane 
(Fig. [5]) that satisfies all three constraints (with IUPAC 
name l,7-diazacyclododecane-4,10-diamine). Lithiation 
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visualizations were generated with PYMOL 3l|. J.E.M. 
thanks Jay Deep Sau for helpful discussions. 
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(1.18) 



FIG. 2. (color online) Top and side views of (a) isolated 
Li and l,7-diazacyclododecane-4,10-diamine, (b) Li bound to 
l,7-diazacyclododecane-4,10-diamine, and (c) Li bound to the 
VN4 defect in diamond. Relative formation energies of Li 
(and Li + in parentheses) are reported in eV, from PBE total 
energies (and PBE-e ionization energy for Li + in (c)). The 
structures for Li + are similar to the neutral structures shown. 



of this molecule should abstract H (as 5H2) to form a 
more stable (by 0.24 eV in our calculations) lithamide. 

Figure [5] depicts lithium in three metastable positions 
relevant to CVD: the isolated atom, bound to a precursor 
molecule, and bound to a B center in diamond. Li + is 
strongly bound to both sites. Li is not bound to the B 
center, but will remain trapped there because of a large 
energy barrier. The Li-N bond length of 2.0 A within the 
precursor is reduced to 1.72 A within the B center. Bond 
strain can be quantified indirectly by comparing total 
energies of the relaxed VN4 cavity to the LiN4 defect 
structure with Li removed, which results in a difference of 
0.27 eV. The reduction of LUNVs ionization energy from 
the precursor to diamond is caused by a destabilization 
of the neutral state rather than any significant change of 
the ionized state. 

Successful CVD synthesis of LiN4 is contingent on the 
existence of growth conditions that preserve the internal 
structure of the LiN4 precursor while still enabling good 
diamond crystal formation, which is an open problem. 

In short, we propose LiN4 as a new donor complex 
in diamond with a predicted activation energy of 0.27 ± 
0.06 eV. Synthesis of LLN4 is likely to require that Li 
be reacted with a preformed VN4 complex, either within 
diamond or a precursor molecule. While further studies 
of Li-N-V defect chemistry in diamond are warranted, 
the present result should serve as sufficient impetus for 
the pursuit of experimental realization of LiN4. 
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Supplemental Material to "Theoretical Design 
of a Shallow Donor in Diamond by Lithium- 
Nitrogen Codoping" 

All DFT calculations are performed using version 5.2 
of the Vienna Ab-initio Simulation Package (vasp) 
A planewave cutoff of 400 eV is used for all calcula- 
tions, along with the manual-recommended [2] projector 
augmented-wave pseudopotentials (or alchemical mix- 
tures thereof Q). Band occupations are set by Gaussian 
smearing with a width of 0.01 eV. Relaxed structures 
have the forces on all atoms reduced below 0.1 eV/A. The 
PBE lattice constant of 3.57 A is used for all calculations 
in diamond at zero pressure. Finite pressure is simulated 
by reducing the lattice constant. A 2 x 2 x 2 Monkhorst- 
Pack grid is used to sample the Brillouin zone (BZ) of 
diamond supercells for PBE total energy calculations, as 
suggested in previous studies Q. All molecules are sim- 
ulated at the T-point of a cubic supercell with a 20 A 
lattice constant. Predictions of molecular structure are 
made by generating many locally stable conformations 
using the MMFF94s force field in AVOGADRO [H] and fur- 
ther relaxing them in VASP. The lowest energy structure 
found by this search is presumed to be the ground state. 

The computational cost of PBE-e calculations is larger 
than PBE because of the Fock exchange step. The ratio 
between costs grows linearly with the number of points 
used to sample the BZ, which limits our PBE-e calcula- 
tions to a single BZ point at a time for large supercells. 
When computing the band structure at a BZ point, the 
density matrix is constructed from only that BZ point 
and iterated to self-consistency. This is a systematic er- 
ror for finite supercells, but it is exponentially suppressed 
with increasing supercell size because all bands are fully 
occupied (half-filled impurity bands are fully occupied in 
one spin channel and empty in the other). 

The addition of Fock exchange in switching from PBE 
to PBE-e has the basic effect of making the impurity 
bands deeper and more sensitive to finite-size effects. 
This is illustrated in Fig. [3] with a plot of the lowest 
branch of the conduction band and the donor impurity 
band of LiN4 in a 6 x 6 x 6 supercell along the L — T — X 
high symmetry path. PBE and PBE-e produce similar 
results, except that the PBE-e impurity band is rigidly 
shifted downwards in energy by 0.68 eV. The average 
depth of the PBE impurity band is smaller than the ex- 
pected shallow impurity depth of 0.2 eV and does not 
change significantly in larger supercells. The PBE-e im- 
purity band is deeper than expected, but it becomes more 
shallow with increasing supercell size and decreasing im- 
purity density. To estimate the isolated impurity limit, 
some kind of extrapolation must be performed. 

We calculate impurity ionization energies A lonlze by 
using a range of diamond supercells from 5 x 5 x 5 to 
8x8x8. Smaller supercells have finite-size effects that 
are too large and complicated to be fit to a simple extrap- 
olation model and larger supercells are computationally 
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FIG. 3. PBE and PBE-e band structures of LiN 4 in a 6 x 6 x 6 
supercell in the spin channel where the donor impurity band 
is occupied. Only the donor impurity band and the lowest 
branch of the conduction band are shown. The conduction 
band is similar at both levels of theory and only one is plotted. 

intractable. We assume a tight binding picture where 
impurity states are localized about each impurity and an 
impurity band manifests from hopping between impu- 
rity states. The hopping energy and impurity bandwidth 
should decrease exponentially with increasing distance 
between neighboring impurities R. All finite-size effects 
in the impurity band should follow a simple exponential 
form cx exp(— R/Rq). Rq is an estimate of the effective 
Bohr radius of the impurity state, but it may be par- 
tially contaminated by finite-size BZ sampling errors. If 
the finite-size dependence is similar and simple enough, 
then some combination of BZ points in the impurity band 
should cancel the finite-size effects. Specifically, we make 
a simple but arbitrary choice to define 

A™ « E C (T) - [(1 - x)E D (T) + xE D (X)] (8a) 
A™ « [(1 - x)E A (T) + xE A (X)} - E V (T), (8b) 

where E c is the lowest branch of the conduction band, 
E v is the highest branch of the valence band, Ed is the 
occupied donor impurity band, and Ea is the unoccupied 
acceptor impurity band. For large supercells, this expres- 
sion will be independent of the choice of x. With this free 
parameter and the results from 4 supercells, we perform 
a least squares fit to minimize the size-dependence of 
A lomze in Eq. (jHJ about an asymptotic ionization energy. 
Additionally, we fit E D (X)-E D (T) and E A (T)-E A (X) 
to the form C exp(— R/Rq) to estimate Rq. 

Extrapolation results are tabulated in Table [TlJ 
(5A lonlzc is the root-mean-square (RMS) deviation be- 
tween the extrapolated A ionizo and Eq. © for the 4 
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TABLE II. Extrapolated donor and acceptor properties. Ro 
is the effective Bohr radius in A. A lonlzc is the ionization en- 
ergy and ^A lonlzc is its RMS variance, both in eV. x is the 
parameter used in the extrapolation formulae in Eq. JSJ. 



Defect 


Ro 






X 


^5 


7.9 


0.31 


0.03 


0.74 


LiN 4 


7.8 


0.27 


0.03 


0.74 


BeN 3 


6.7 


0.39 


0.04 


0.69 


BN 2 


5.5 


0.50 


0.03 


0.62 


P 


5.2 


0.54 


0.02 


0.61 


N 


4.5 


0.86 


0.04 


1.04 


B 


7.9 


0.31 


0.03 


1.05 


cf 


7.0 


0.30 


0.01 


0.95 



calculated supercells. The extrapolation procedure suc- 
ceeds in reducing the variance to a value significantly 
below the bandwidth of the impurity bands in the calcu- 
lated supercells. The effective Bohr radii extracted from 
extrapolation increase with decreasing ionization energy, 



as expected from a simple hydrogenic model of impuri- 
ties. The hydrogenic impurity model gives a radius of 
6.3 A for shallow donors and 2.8 A for shallow acceptors, 
which is smaller than the extrapolated values and espe- 
cially so for the acceptors. The cause of this deviation is 
unclear, but the effect on the ionization energies seems 
minor since boron is within 0.06 eV of its experimental 
activation energy of 0.37 eV. 
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